Homotopy analysis of 1D unsteady, nonlinear groundwater flow through porous media

SONG, H. and TAO, L., 2007. Homotopy analysis of 1D unsteady, nonlinear groundwater flow through porous media. Journal of Coastal Research, SI 50 (Proceedings of the 9th International Coastal Symposium), pg – pg. Gold Coast, Australia, ISBN 1-891276-54-9 In this paper, the 1D unsteady, nonlinear groundwater flow through porous media, corresponding to flood in an aquifer between two reservoirs, is studied by mass conservation equation and Forchheimer equation instead of Darcy's law. The coupling nonlinear equations are solved by homotopy analysis method (HAM), an analytic, totally explicit mathematic method. The method uses a mapping technique to transfer the original nonlinear differential equations to a number of linear differential equations, which does not depend on any small parameters and is convenient to control the convergence region. Comparisons between the present HAM solution and the numerical results demonstrate the validity of the HAM solution. It is further revealed the strong nonlinear effects in the HAM solution at the transitional stage. ADDITIONAL INDEX WORDS: Homotopy analysis method, Forchheimer equation, porous media INTRODUCTION Groundwater flow through porous media is traditionally described by Darcy’s law and it is normally valid for low Reynolds pore-scale numbers. For moderate and high velocity flow, however, Forchheimer equation should be applied due to nonlinear effects. STARK (1972) numerically solved the NavierStokes laminar flow equations and tested the relations of Darcy’s law, Forchheimer equation and others; INNOCENTINI et al. (1999) compared Darcy’s law and Forchheimer equation and recommended highly the latter in order to take into account permeability; NIELD (2000) discussed the inertial effects on viscous dissipation for the case of Darcy, Forchheimer and Brinkman models. Forchheimer equation could be derived in different approaches (e.g. AHMED and SUNADA, 1969; HASSANIZADEH and GRAY, 1987; WHITAKER, 1996) and has been proved in theoretical and experimental way (MACDONALD et al., 1979; THAUVIN and MOHANTY, 1998). Extensive studies on the parameters in Forchheimer equation have been carried out. COULAUD et al. (1988) introduced a nonlinear term into Darcy’s equation and solved it numerically. In his approach, the hydrodynamic constants in the Forchheimer equation were expressed by the expression of porosity and geometrical data. WANG and LIU (2004) investigated the scaling relations for the fluid permeability and the inertial parameter in the Forchheimer equation, by solving the Navier-Stokes equation for flow in a two-dimensional percolation porous media. Although the application of the Forchheimer flow is very useful and practical, very limited attempts on the analytical approach have been reported in the literature. MOUTSOPOULOS and TSIHRINTZIS (2005) solved the Forchheimer flow through porous media in 1D form by perturbation method, dividing the problem into two stages and solving them by two sets of equations. For numerical method, GREENLY and JOY (1996) used onedimensional finite element method and Forchheimer equation to investigate the groundwater flow through a valley fill. EWING et al. (1999) used finite difference, Galerkin finite element and mixed finite element techniques to investigate Forchheimer flow in a hydrocarbon reservoir. KIM and PARK (1999) and PARK (2005) used mixed finite element method to analyse the flow of a singlephase fluid in a porous medium governed by Forchheimer equation. Recently, a new mathematical technique, namely homotopy analysis method (HAM) has been applied to nonlinear fluid dynamics problems (LIAO, 1995, 2004). The approach does not depend on small or large parameters and is easy to adjust the convergence region and rate of approximation series. In this paper, homotopy analysis method is applied to solve the 1D unsteady, nonlinear groundwater flow through porous media, corresponding to a flood in a long aquifer between two reservoirs, or to a flow in a laboratory column restrained by two external tanks. The coupling equations are transformed by similarity law and a global solution, which is analytical, totally explicit is obtained. The piezometric head from the present HAM solution for nonlinear flow agrees well with numerical results and the previous perturbation solutions. THEORETICAL CONSIDERATION Consider a one-dimensional (1D) flow in a confined porous medium shown in Figure 1. Before t=0, the piezometric head is a constant hI. After t=0, the piezometric head at left end raises Δh, while the piezometric head at the far right end remains unchanged. Journal of Coastal Research, Special Issue 50, 2007 Homotopy analysis of the 1D non-steady, nonlinear groundwater flow through porous media The movement of the flow satisfies the equation of mass conservation as following: ( ) h S Bq t ∂ +∇ = ∂ r R , (1) where S is the storage coefficient; h is the piezometric head; t is the time; ∇ is the 2D Nabla operator; B is the thickness of the aquifer; is the velocity; and R is the external sink-source term, which is assigned zero in this particular problem. q Assuming that the properties of the aquifer are homogeneous, the Forchheimer or Forchheimer-Dupuit equation is: h aq bq q −∇ = + | | r r r , (2) where a and b are coefficients. Equations (1) and (2) can be expressed in 1D form as 0 h q S B t x ∂ ∂ + = ∂ ∂ , (3) 2 h aq bq x ∂ − = + ∂ . (4) The initial condition is: at 0 I h h t = = . (5) The boundary conditions are: for 0 at 0 I h h h x t = + Δ = > , (6) for at any I h h x t = = +∞ . (7) Introducing I h h h h − = Δ % , (8) Equations (1) and (2) are transformed as: 0 S h q h B t x ∂ ∂ Δ ⋅ + = ∂ ∂ % , (9) 2 h h aq bq x ∂ −Δ = + ∂ % . (10) From Equation (10), the velocity can be expressed as 2 4 2 h x a a b h q b ∂ ∂ − + − Δ = % . (11) If 2 4 x b h a ∂ − Δ > % , the inertial term is dominant, otherwise the Darcy (viscous) term is dominant. Substituting Equation (11) into Equation (9), we have: 2 2 1 2 2 h h h C x t x ⎛ ⎞ ∂ ∂ ∂ C / = + ⎜ ⎟ ∂ ∂ ∂ ⎝ ⎠ % % % , (12) subject to the initial and boundary conditions: 0 for h = % 0 t = 0 x , (13) 1 for h = = % , (14) 0 for h x = = +∞ % , (15) where ( ) 1 4 S B C b h = − Δ , ( ) 2 2 S B = C a . Using the similarity transformations: 1 2 2 ( ) ( ) h x t t f x t ξ ξ / , = , = % / , (16) Equation (12) becomes 2 1 2 1 2 ( ) 2 ( ) 16 2 ( ) ( ) 2 ( ) f f C f C f f ξ ξ ξ ξ ξ ξ ξ ξ / ′ ′′ ⎛ ⎞ + ′ ⋅ = + ⎜ ⎟ ′ − ⎝ ⎠ , (17) with boundary conditions: 3 (0) f C = , (18) ( ) 0 f +∞ = , (19) where 1 2 3 1 C t / = / , which is a constant for a given time t. Alternatively, Equation (17) can be expressed as,